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Abstract 

We have developed a general technique to study the dynamics of the quantum adiabatic evolution 
algorithm applied to random combinatorial optimization problems in the asymptotic limit of large 
problem size n. We use as an example the NP-complete Number Partitioning problem and map 
the algorithm dynamics to that of an auxilary quantum spin glass system with the slowly varying 
Hamiltonian. We use a Green function method to obtain the adiabatic eigenstates and the minimum 
excitation gap, g m { n = 0(n2~ n / 2 ), corresponding to the exponential complexity of the algorithm 
for Number Partitioning. The key element of the analysis is the conditional energy distribution 
computed for the set of all spin configurations generated from a given (ancestor) configuration by 
simulteneous tipping of a fixed number of spins. For the problem in question this distribution is 
shown to depend on the ancestor spin configuration only via a certain parameter related to the 
energy of the configuration. As the result, the algorithm dynamics can be described in terms of 
one-dimenssional quantum diffusion in the energy space. This effect provides a general limitation 
on the power of a quantum adiabatic computation in random optimization problems. Analytical 
results are in agreement with the numerical simulation of the algorithm. 
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I. INTRODUCTION 



Since the discovery by Shor |T[ nearly a decade ago of a quantum algorithm for effi- 
cient integer factorization there has been a rapidly growing interest in the development of 
new quantum algorithms capable of solving computational problems that are practically 
intractable on classical computers. Perhaps the most notable example is that of a combina- 
torial optimization problem (COP). In the simplest case the task in COP is to minimize the 
cost function ("energy") E z defined on a set of 2 n binary strings z = {zi, . . . , z n } Zj = 0,1, 
each containing n bits. In quantum computation this cost function corresponds to a Hamil- 
tonian Hp 

H p = J2e z \z)(z\ (1) 

z 

|z) = |^i)i <g) (2:2)2 <£>•••<£> \z n ) n - 

where Zj = 0, 1 and the summation is over 2 n states |z) forming the computational basis 
of a quantum computer with n qubits. State \zj)j of the j-th qubit is an eigenstate of the 
Pauli matrix a z with eigenvalue Sj — 1 — 2zj (Sj = ±1). It is clear from the above that the 
ground state of Hp encodes the solution to the COP with cost function E z . 

COPs have a direct analogy in physics, related to finding ground states of classical spin 
glass models. In the example above bits Zj correspond to Ising spins Sj. The connection 
between the properties of frustrated disordered systems and the structure of the solution 
space of complex COPs has been noted first by Fu and Anderson [Q. It has been recognized 
H that many of the spin glass models are in almost one-to-one correspondence with a 
number of COPs from theoretical computer science that form the so-called NP-complete 
class @] . This class contains hundreds of the most common computationally hard problems 
encountered in practice, such as constraint satisfaction, traveling salesmen, and integer 
programming. NP-complete problems are characterized in the worst cases by exponential 
scaling of the running time or memory requirements with the problem size n. A special 
property of the class is that any NP-complete problem can be converted into any other NP- 
complete problem in polynomial time on a classical computer; therefore, it is sufficient to 
find a deterministic algorithm that can be guaranteed to solve all instances of just one of the 
NP-complete problems within a polynomial time bound. It is widely believed, however, that 
such an algorithm does not exist on a classical computer; whether it exists on a quantum 
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computer is one of the central open questions. Ultimately, one can expect that the behavior 
of new quantum algorithms for COPs and their complexity will be closely related to the 
properties of quantum spin glasses. 

Recently, Farhi and co-workers suggested a new quantum algorithm for solving combina- 
torial optimization problems which is based on the properties of quantum adiabatic evolution 
||. Running of the algorithm for several NP-complete problems has been simulated on a 
classical computer using a large number of randomly generated problem instances that are 
believed to be computationally hard for classical algorithms || |7], ||, ||. Results of these 
numerical simulations for relatively small size of the problem instances ( n < 20) suggest 
a quadratic scaling law of the run time of the quantum adiabatic algorithm with n. Fur- 



thermore, it was shown in [10] that the previous query complexity argument that led to 



the exponential lower bound for unstructured search [11] cannot be used to rule out the 
polynomial time solution of NP-complete Satisfiability problem by the quantum adiabatic 
algorithm. 



In |L0|, [12], [13], [14], [T5|] special symmetric cases of COP were considered where symmetry 
of the problem allowed the authors to describe the true asymptotic behavior (n — > oo) of the 
algorithm. In certain examples considered in 0, [U| the quantum adiabatic algorithm finds 
the solution in time polynomial in n while simulated annealing requires exponential time. 
This effect occurs due to the special connectivity properties of the optimization problems that 
lead to the relatively large matrix elements for the spin tunneling in transverse magnetic 
field between different valleys during the quantum adiabatic algorithm. In the examples 
considered in |13|] the tunneling matrix element scales polynomially with n. On the other 
hand, in simulated annealing different valleys are connected via classical activation processes 
for spins with probabilities that scale exponentially with n. It was also shown for certain 
simplified examples |L4|, [L5| , that quantum adiabatic algorithm can be modified to completely 
suppress the tunneling barriers even if the corresponding classical cost function has local 
minima well separated in the space of spin configurations. 

However, so far there are no study on the true asymptotic behavior of the algorithm 
for the general case of randomly generated hard instances of NP-complete problems. Also 
there are no analysis of the limitations of the quantum adiabatic computation arising from 
the intrinsic properties of disorder and frustration in this problems. Such analysis is of the 
central interest in this paper. 
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In Sec. |T| we introduce the random Number Partitioning problem and describes condi- 
tional cost distributions (neighborhood properties) in this problem. In Sec. [ITT] we describe 
the concept of quantum adiabatic computation applied to combinatorial optimization prob- 
lems and introduce a Green function method for the analysis of the minimum gap. In Sec. [TV] 
we describe the effect of quantum diffusion in the algorithm dynamics, derive the scaling for 
the minimum gap and the complexity of the algorithm for the random Number Partitioning 
problem. We also obtain the scaling of the minimum gap numerically from the form of the 
cumulative density of the adiabatic eigenvalues at the avoided-crossing point. In Sec. [V] we 
discuss the results of the simulations of the time-dependent Schrodinger equation to sim- 
ulate quantum adiabatic computation for Number Partitioning and obtain its complexity 
numerically. 



II. NUMBER PARTITIONING PROBLEM 



Number Partitioning Problem (NPP) is one of the six basic NP-complete problems that 
are at the heart of the theory of NP-completeness It can be formulated as a combinatorial 
optimization problem: Given a sequence of positive numbers {a\, . . . , a n } find a partition, 
i.e. two disjoint subsets A and A', such that the residue 



E 



Oj eA a-j eA' 

is minimized. In NPP we search for the bit strings z = {z\, . . . , z n } (or corresponding Ising 
spin configurations S = {Si, . . . , S n }) that minimize the energy or cost function E x 



n 



^ajSj, Sj 

3=1 



1 - 2z 



■J- 



(3) 



where Sj = 1 [zj = 0) if aj e A and Sj = —1 (zj = 1) if ctj G A'. The partition S 
with minimum residue can also be viewed as the ground state of the Ising spin glass, — 0|, 
corresponding to the Mattis-like antiferromagnetic coupling, J^- = — ajOj. 



NPP has many practical applications including multiprocessor scheduling []T6j| , cryptogra- 



phy [ IT | , and others. The best deterministic heuristical algorithm for NPP, the differencing 
method of Karmakar and Karp |Tj|, can find with high probability solutions whose energies 
are of the order \/n aXogn for some a > 0. The interest in NPP also stems from the re- 
markable failure of a standard simulated annealing algorithm for the energy function (^]) to 



find good solutions, as compared with the solutions found by deterministic heuristics |L9 . 
The apparent reason for this failure is due to the existence of order 2 n local minima whose 
energies are of the order of 1/n [pO] which undermines the usual strategy of exploring the 
space of the spin configurations S through single spin flips. 

The computational complexity of random instances of NPP depends on the number of bits 
b needed to encode the numbers aj. In what follows we will analyze NPP with independent, 
identically distributed (i.i.d.) random b-bit numbers aj. Numerical simulations show |2l|, [22|. 
[26| that solution time grows exponentially with n for n <C b then decreases steeply for n ~>b 
(phenomenon of "peaking") and eventually grows polynomially for n ^> b. The transition 
from the "hard" to computationally "easy" phases at n ~ b has features somewhat similar 
to phase transitions in physical systems |23[ . The detailed theory of the phase transition in 



NPP was given in Refs. [24, |25|]. If one keeps the parameter £ = b/n fixed and lets n — ► oo 
then instances of NPP corresponding to £ > 1 will have no perfect partitions with high 
probability. On the other hand for £ < 1 number of perfect partitions will grow exponentially 
with n. Transitions of this kind were observed in various NP-complete problems |p8[| . In 
what follows we will focus on the computationally hard regime £ 3> 1. 

A. Distribution of signed partition residues 

The values of individual energies are random and depend on the particular instance of 
NPP (i.e., the set of numbers aj). However on a coarse-grained scale (i.e. after averaging 
over individual energy separations) the form of the typical energy distribution is described 
by some universal function for randomly generated problem instances. We introduce for a 
given set of randomly sampled numbers aj a coarse-grained distribution function of signed 
partition residues Q z ([|) 

P(ty = 2- n — / drj ( 4 ) 

Jn-AQ/2 ze {o,l}« 

Here S(x) is the Dirac delta-function; the sum is over 2 n bit-strings z and 2~ n is a normal- 
ization factor. In (M we average over an interval AQ of the partition residues whose size is 
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(5) 



chosen self-consistently, AQ ^> 2 n /P(fi). Using (^) we can rewrite ([|) in the form 

i pco /AO \ 

P(O) = — jf rf S C J /(a) cos(^), 

n 

I(s) = JJcos(ajs), C(a;) = sin(x)/x. 
i=i 

Here ((x) is a window function that imposes a cut-off in the integral @ at s ~ 2/ All For 
large n this integral can be evaluated using the steepest descent method. In the following 
we shall assume that the 6-bit numbers aj are distributed inside of the unit interval [0, 1] 
and are integer multiples of 2~ b , the smallest number that can be represented with available 
number of bits b. We note that for large n the function I(s) has sharp maxima (minima) 
with width ~ n^ 1 ^ 2 at the points Sk = k7i2 b , k = 0,1,...; |/(s/c)| = 1. Only one saddle 
point at s = contributes to the integral in (f|) due to coarse-graining of the distribution (|]). 
Indeed, it will be seen below that the window size 2/AQ can be chosen to obey the conditions 
1 < n l / 2 /AVL < 2 n . Therefore in the case of high-precision numbers, &>b, saddle-points 
Sk with k > lie far outside the window and their contributions can be neglected (see also 
Appendix On the other hand the window function C(x) can be replaced by unity while 
computing the contribution from the saddle-point at s = 0. Finally we obtain for \Q\ n 

( C f. a) 



n 

a 2 (0) = -^a3 (£<n). (6) 

71 3=1 

The coarse-grained distribution -P(fi) depends on the set of a/s through a single self- 
averaging quantity <r(0) (cf. p3|). 

One can also introduce the distribution P(E) of cost values (energies) E z = |f2 z |. Due to 
the obvious symmetry of the NPP, the cost function E z in (|3|) does not change after flipping 
signs of all spins, Sj — > —Sj. Therefore 

P(E) = 1/2P(±E). (7) 

We emphasize that, according to Eq. @ for a typical set of high-precision numbers a,j the 
energy spectrum in NPP is quasi-continuous, and there are only two scales present in the 
distribution P{E)\ one is a "microscopic" scale given by the characteristic separation of the 



individual partition energies, E min , and another is given by the mean partition energy (E) 
(or the distribution width (E 2 ) 1 ^ 2 ) 

E min ~ a(0) n 1 ' 2 2- n , (E 2 ) = |(£> 2 = na 2 (0). (8) 

This justifies the choice for AQ above that corresponds to coarse-graining over many indi- 
vidual energy level separations. 

We note that the distribution P(fi) (|6]) is Gaussian for E <C n and can be understood 
in terms of a random walk with coordinate Q using Eq. ([3]). The walk begins at the origin, 
Q = 0, and makes a total of n steps. At the j-th step Q moves to the right or to the left 
by "distance" 2aj if Sj = 1 or Sj = — 1, respectively. In the asymptotic limit of large n the 
result @ corresponds to equal probabilities of right and left moves and the distribution of 
step lengths coinciding with that of the set of numbers {2aj}. 

Finally, the energy distribution function P(E) of the form ©,(0) was previously obtained 
by Mertens using explicit averaging over the random instances of NPP. He also computed 
the partition function Z(T) for a given instance of NPP at a small finite temperature T using 
the steep est- descent method and summation over the saddle-points Sk = kir 2 b similar to our 



discussion above [23J (in his analysis T played a role similar to our regularization factor 
AQ in (|),(|)). 



We emphasize however, that the approach in Ref. [23] based on Z(T) is necessarily 
restricted to the analysis of the "static" properties of NPP at E ~ 2~ n , i.e., the phase 
transition in the number of perfect partitions |23j when the control parameter £ = n/b 
crosses a critical value. On the other hand distribution P(O) @) introduces at finite energies, 
as well as the conditional distribution introduced in the next section also allow us to directly 
study the intrinsic dynamical properties of the problem in question such as the dynamics of 
its quantum optimization algorithms. 

B. Conditional distribution of signed partition residues 

Consider the set of bit-strings z' obtained from a given string z by flipping r bits. The 
conditional distribution of the partition residues Q z > (|3|) in the r-neighborhood of z can be 
characterized by its moments: 

<° fc > = Cf) E (^') fc W, Z ); k = l,2,... (9) 



Here S mi i is a Kronecker delta and function D(z, z') computes the number of bits that take 
different values in the bit-strings z and z'. It is the so-called Hamming distance between the 
two strings 

n 

D(z,z')= J>,--^|. (10) 

3=1 

The Hamming distance r = D(z, z') between the bit-strings is directly related to the overlap 
factor q between the corresponding spin configurations often used in the theory of spin 
glasses |, fZjJ: 

1 n 2 
q=-T,S j S' j = l--D(*,*). (11) 

3=1 

(in what follows we shall use both quantities r and q). For A; = 1, 2 in @ one obtains after 
straightforward calculation the first and second moments of the conditional distribution 

(fi)=gfi., (12) 
<fi 2 > - (fi> 2 = ^ 2 (g) f 1 + -±-) f 1 - I -g-V (13) 



n — 1/ \ n (P 2 ) 
a(g) = cx(O) (1 - q 2 fl\ q= 1 - ^, (14) 

where cr(0) and (P 2 ) are given in (|6]) and (|j), respectively. 

The conditional distribution of Q z > can also be defined in a way similar to (f|) 

Try/ ^ E *fa " «-0 <W») (15) 

W ^« JSl'-Ml'/2 Z 'g {0 ,l}n 

where averaging is over the small interval AQ' that, however, includes many individual 
values of Q, z > for a given r. It is clear from (|i~2f ) , (|i~3l) that the first two moments of P r}Z (Q') 
depend on z only via the value of Q z . This does not hold true, however, for the higher- 
order moments that depend on other functions of z as well. For example, (Q 3 ) involves the 
quantity Y^j=x a f(l ~~ 2zj), etc. 

Our main observation is that in the asymptotic limit of large n the conditional distri- 
bution P r z (fl') is well-described by the first two moments (|T2|) , ({Op . Then, according to 
the discussion above, its dependence on z is only via Q z . The detailed study of the higher 
moments (§) will be done elsewhere. Here we use the following intuitive approach relevant 
for analysis of the computational complexity of the quantum adiabatic algorithm for the 
NPP. We average P rz (f2 / ) over the strings z with residues Q z inside a small interval AQ 
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(containing, however, many levels Q z ). After such averaging the result, P r (Q'\Q), can be 
written in the form 

Pr(P\n) = ^ffi, (16) 
i pQ+AQ/2 

p r (n',n) = 2- n — dr, J2 s( v -n z )p r , z (n'), (17) 

where -P(O) is given in (|6]). We note that Eq. ( |T6| ) formally coincides with the Bayesian 
rule expressing the conditional distribution function P r (Q'\Q) through the 2-point (joint) 
distribution function P r (Q',Q) and the single-point distribution P(O) (^). 

The explicit form of P r (Q\Q) is derived in Appendix [B] in a manner similar to the deriva- 
tion of P{E) in Sec. [II A| . The results are presented in Eqs. (|B7|) and ( p8|) . They show that 
P r (Q',Q) in the limit n 3> 1 is indeed well described by its first two moments that corre- 
sponds precisely to the expressions given in Eqs. (|12|)'(0) above. From this we conclude 
that 

P rtZ (Q') = P r (n'\Q z ). (18) 

In the case r = 1 there are n strings z' at a Hamming distance 1 from the string z. 
Partition energies corresponding to these strings equal |flg — 2a,jSj\, 1 < j < n (cf. (§)). 
After the coarse-graining over the energy scale 0(l/n) in the range, \Q'\ <C n, the 
conditional distribution P rz is a step function in the interval Q z — Q' G [—2, 2]. For r = n — 1 
one has the same form of the distribution but for fl z + Q'. Both results correspond to 
nearly equal distribution of spins between between ±1 values. Then in the range of energies 
1^1 j l^z| ^ 1 one has: 

P r , z (n')~ P r = l/2 + , r = l,n-l (n>l) (19) 

For r,n — r ^> 1 distribution P rz (f2') has a Gaussian form with a broad maximum at 
Q' = qQ z (cf. Eqs. ([T2|) , (p~3|) , 7|) ) . Near the maximum we have: 

P, Z (Q') « P r = 1 \Q'\, \Q Z \ « r^Vfe). (20) 

We studied the conditional distribution in NPP numerically as well (see Fig.|TJ and Sec.[B|). 
The results are in good agreement with theory even for modest values of n < 30. 

The characteristic spacing between the values of the partition residues in the subset of 
strings z' with D(z',z) = r is l/(P r (")) f° r n °t too large E z ,E z i (see above). This spacing 



1.5 



s 



1 



I 



i 



V 




I 



4 



26 



r 



FIG. 1: Plots of the (scaled) conditional distribution |jj) s = a(0)(2vrn) 1 / 2 (Afi)" 1 Jjf drj P r , z (v) 
vs r are shown with points. We use coarse-graining window AS7=0.3. Different plots correspond to 
29 randomly selected bit-strings z with energies |f2 z | G [0,0.3] for one randomly generated instance 
of NPP with n = 30 and 6 = 35. For r,n-r> 1 the values of s corresponding to different strings 
are visually indistinguishable from each other. Dashed line is a plot of a(0)/a(q) vs r given in 
( |l3| ) (q = 1 — 2r/n). Insert: plots of the integrated quantity given in (fBg|), Q = ^ J drj P r , z {v) 
vs x = Q/(a(q)\^2n), for different values of r = 2, . . . , ra/2 and randomly selected bit-string z with 
energy |f2 z | close to 0. All plots correspond to the same instance of NPP as the main figure. Plots 
for different values of r are visually indistinguishable from each other and from the theoretical 
curve given in ( [BIO] ). 

decreases exponentially with the magnitude of the string overlap factor, |g| = \(n — 2r)/n\. 
The hierarchy of the subsets corresponding to different values of |g| form a specific structure 
of NPP. We note that the distribution of partition residues within the hierarchy is nearly 
independent of the ancestor string z in a broad range of energies E' < n 1 ' 2 where P r>z (E') ss 
P r . One can see that the magnitude of the overlap factor q between two strings with energies 
within a given interval [0, E] is limited by some typical value q satisfying the following 
equation: 



The smaller E is, the smaller \q\ is: strings that are close in energy are far away in the 
configuration space. This property gives rise to an exponentially large number of local 
minima for small values of E z that are far apart in the configuration space. For example, 
strings with E z ~ -E m in typically correspond to \q\ = 0(l/n), they can be obtained from 




r 

= 1-2-. 

n 



(21) 
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each other only by simultaneously flipping clusters with ~ n/2 spins. 

Eq. (^) describes the dynamics of a local search heuristic (e.g., simulated annealing). 
It shows that the average cost value E during the search decreases no faster than 0(1/M) 
where M — O ( ( n J J is the number of generated configurations. This result coincides with 
that obtained in |2Dj using a different approach. It says that any classical local search 
heuristic for NPP cannot be faster than random search. Indeed, during local search the 
information about the "current" string z with E z < 1 is being lost, on average, after one 
spin flip (cf. Eqs. (|19|) , (|20|) ) . We show below that precisely this property of NPP also leads 
to the complexity of the quantum adiabatic algorithm corresponding to that of a quantum 
random search. 

We note that one can trivially break the symmetry of NPP mentioned above by introduc- 
ing an extra number and placing it, say, in the subset A. In this case different partition 
energies will still be encoded by spin configurations S = {Si, . . . , S n } (or corresponding bit- 
strings z) with Qs = ^0 + J2j=i Sj a j an d E z = \Qs\ (cf- |3|)- We shall adopt this approach in 
the analysis of the performance of the quantum adiabatic algorithm for NPP given below. 



III. QUANTUM ADIABATIC EVOLUTION ALGORITHM 

In the quantum adiabatic algorithm || one specifies the time-dependent Hamiltonian 
H(t) = H(t/T) 

H(t) = (1-t)V + tH p , (22) 

where r = t/T is dimensionless "time". This Hamiltonian guides the quantum evolution 
of the state vector \il>(t)) according to the Schrodinger equation id\ip{t))dt = H(t)\ip(t)) 
from t = to t = T, the run time of the algorithm (we let h = 1). Hp is the "problem" 
Hamiltonian given in ([TJ) . V is a "driver" Hamiltonian, that is designed to cause transitions 
between the eigenstates of Hp. In this algorithm one prepares the initial state of the system 
ip(0) to be the ground state of -f^(O) = V. In the simplest case 

n 

V = ~Y,< hA(0)> = 2-/ 2 5>>, (23) 

j=l z 

where a J x is a Pauli matrix for j-th qubit. Consider instantaneous eigenstates \<P v (t)) of 
H(t) with energies A r? (r) arranged in nondecreasing order at any value of r 6 (0, 1) 

H\<j )v } = \ ri \ ( f )v ), r / = 0,l,...,2"-l. (24) 
11 



Provided the value of T is large enough and there is a finite gap for all t G (0, T) between 
the ground and excited state energies, g(r) = Ai(r) — Ao(t) > 0, quantum evolution is 
adiabatic and the state of the system \i[)(t)) stays close to an instantaneous ground state, 
|0o(t/T)) (up to a phase factor). Because H(T) = Hp the final state \ip(T)) is close to the 
ground state \4>o(t = 1)) of the problem Hamiltonian. Therefore a measurement performed 
on the quantum computer at i = T (r = 1) will find one of the solutions of COP with large 
probability. 

There is a broad class of COPs from theoretical Computer Science where the number of 
distinct values of a cost function scales polynomially in the size of an input n. An example 
is the Satisfiability problem in which the cost E z of a given string z equals the number 
of constrains violated by the string. For those problems, the spectrum of H(t), at the 
beginning (r ~ 0) and at the end (r pa 1) of the algorithm, consists of a polynomial number 
of well-separated energy levels. Quantum transitions away from the adiabatic ground state 
occur most likely near the avoided-crossing points t pa r* where the energy gap g{r) reaches 
its minima ||. Near the avoided-crossing points, the spectrum of H(r) is quasi-continuous, 
with the separation between individuals eigenvalues scaled down with n. The probability of 
a quantum transition, 1 — \ (ip{t)\(f) (t /T))\^ =T , is small provided that 

T> 5 ' 9mm = mm [Ai(r) - A (r)J , (25) 

Q ■ 0<t<1 
i<min — — 

(H T = dH/dr). The fraction in (p5|) gives an estimate for the required runtime of the 
algorithm and the task is to find its asymptotic behavior in the limit of large n> 1. The 
numerator in (p5|) is less than the largest eigenvalue of H T = Hp — V, typically polynomial 
in n 0. However, g m j n can scale down exponentially with n and in such cases the runtime 
of the quantum adiabatic algorithm will grow exponentially with the size of COP. 



A. Implementation of QAA for NPP 

As suggested in || the quantum adiabatic algorithm can be recast within the conventional 
quantum computing paradigm using the technique introduced by Lloyd f30| . Continuous- 



time quantum evolution can be approximated by a time-ordered product of unitary op- 
erators, e~ l ^ l ~ Tk ^ v& e~ lTkHpS , corresponding to small time intervals (tk,tk + 5). Operator 
e -j 0--r k )V6 typically corresponds to a sequence of 1- or 2-qubit gates (cf. (^3|)). Operator 
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e -iT k H P 8 j g diagonal j n ^ ne computational basis |z) and corresponds to phase rotations by 
angles E z 5. Since in the case n -C b, the average separation between the neighboring values 
of E z is 1/P(E) = 0(2~ n ), the quantum device would need to support a very high precision 
in its physical parameters (like external fields, etc.) to control small 0(2~ n ) differences in 
phases. Since this precision scales with n exponentially it would strongly restrict the size 
of an instance of NPP that could be solved on such a quantum computer. This technical 
restriction is generic for COPs that involve a quasi-continuous spectrum of cost-function 
values. Among the other examples are many Ising spin glass models in physics (e.g., the 
Sherrington-Kirkpatrick model ||). To avoid this restriction we introduce a new oracle-type 
cost function £ z that returns a set of values 



£ z = c(tt z ), c(x) — > {e , ei, . . . , em) (e*+i > e k ), (26) 

that can be stored using a relatively small number of bits O(logn). For example, we can 
divide an interval of partition energies (0,5), B = X^=o a i ^ ms wnose sizes grow 
exponentially with the energy. Then the new cost will take one value per bin 

c(x) = £ k = —M + k for u k < \x\ < u k+ i, 

oo k = (2 k - 1) A, k = 0, . . . , M. (27) 

The last bin is % < \Q Z \ < B where we have £ z = Em = 0. The value of the cutoff ujm < B 
is discussed below. In this example the Hilbert space of 2 n states |z) is divided into M + 1 
subspaces C k , each determined by Eq. ( p7|) for a given k 

M 

ffj> = 5>*X>X*l- (28) 
fc=o ze£ fe 

Note that subspace Cq contains the solution(s) to NPP. Dimension d of Cq is controlled 



by the value of A in (|27 ) which is another control parameter of the algorithm. We set 
A = 2~ n K/P(0) where the integer K w do ^> 1 is independent of n and determines how 
many times on average one needs to repeat the quantum algorithm in order to obtain the 
solution to NPP with probability close to 1. 

Operator Hp projects any state onto the states with partition residues in the range 
< |O z | < ujm- If we choose 

l<u M <^(E), (29) 
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then the distribution function (|6]) is nearly uniform for \Q Z \ < ujm- Therefore the dimensions 
of the subspaces grow exponentially with k: dk = do 2 k for k < M. This simplification 
would not affect the complexity of a quantum algorithm that spends most of its time in "an- 
nealing" the system to much smaller partition residues, ujm 3> |^ z | ~ B m m = Oir^l 2 2~ n ). 

We note that the new discrete-valued cost function defined in ( 27 ) is non-local. Unlike 
problems such as Satisfiability, it cannot be represented by a sum of terms each involving a 
small number of bits. To implement a unitary operator e~ tTkHpS with Hp given in one 
needs to implement the following classical function on a quantum computer 



£ z = Q(uj m - |^z|) 



log 2 



Here [x] denotes the integer part of a number x; Q(x) is the theta-function (Q(x) = 1 for 
x > and Q(x) = for x < 0). The implementation of (|30|) with quantum circuits involves, 
among other things, the addition of n numbers together with their signs to compute Q z , and 
taking the discrete logarithm of a 6-bit number with respect to base 2. These operations 
can be performed using a number of quantum gates that is only polynomial in n and b (cf. 
|]J for the implementation of the discrete logarithm). 

Since the implementation of a cost function (|25|) , fl5H|) does not add an exponential over- 
head to the complexity of QAA the feasibility of this algorithm for NPP depends on the 
scaling of the minimum gap g m i n with n. 



B. Stationary Schrodinger equation for adiabatic eigenstates 

We now solve the stationary Schrodinger equation ( f24"|) and obtain the minimum gap g min 
( p5D in the asymptotic limit n — > oo. To proceed we need to introduce a new basis of states 
l x ) = \ x i)i ® ^2)2 <8> • • • <8> \x n ) n where state \xj)j is an eigenstate of the Pauli matrix a x for 
the j-th qubit with eigenvalue 1 — 2xj — ±1. Driver Hamiltonian V can be written in the 
following form: 

n 

m=0 x\-\ \-x n =m 

For a particular case given in Eq. (^) we have V m = 2m — n. Matrix elements of T m in a 
basis of states |z) depend only on the Hamming distance D(z,z') between the strings z and 
z' 

(z|x-|z') = rs {sy) , (32) 

14 



C = 2- EE (% Q(-l) P A m „ +P . (33) 



g=0 p=0 

We now rewrite Eq. fl2~4|) in the form 

|0) = -^_ F P |0), a = a(r) = l-r, (34) 

A — OJK 

(we drop the subscript t] indicating the number of a quantum state and also the argument r 
in <p and A). From (|2?|)-(|3~4]) we obtain the equation for the amplitudes <p z = (z|0) in terms 
of the coefficients I™ 

[1 - tG c(fi z )] Z = — IrV Gd^) z / c(fi z ), (35) 

A — a Vo 77 



$ = E C (Q 



z' / fz' j 

z' 



rm 

G r = G r (X) = V x r T/ , < r < n. 



Here we separated out a "symmetric" term oc 2 n $ corresponding to the coupling between 
the states |z) via the projection operator X° 



IV. MINIMUM GAP ANALYSIS 

A. Coarse-graining of the transition matrix 



We now make a key observation that Z in fl35[) can be determined based on the properties 



of the conditional distribution P r ^(E) (]15|) and the form of the Green function G r (X). We 
sum the Green function Gd( z ,z') over all possible transitions from a given state z' to states 
z' 7^ z with energy e^. For not too large partition residues of the initial and final states we 
obtain 

E G DM {\) « F k {\) + / z ,, fc (A) (36) 

zGjC^ , z^z' 

Ft(A) = ^,. W = ^*_^_Q Gr(A) (37) 

|fi z ,|, |O z | « (£>, /i=^- (38) 

Function <j(g) above is defined in flT4]) and / z ',fc(A) is a small correction described below. In 
function s(A) we replaced summation over the integer values of r by an integral. It can be 
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evaluated using the explicit form of G r (X) that decays rapidly with r. In what follows we 
will be interested in the region |A — aVo\ <C 1 where 

^ (m+r) 

m=l 

(7 is Euler's constant) and s(A) ~ — ln2/(2a). We note that 



(v -1 n-r Q—n I n \ 
n\ y UW _ 2 -"(lnr + 7 ). (39) 
r / m 



-2aG r (\) « - 



(n/2 - r) 



-1 



n/2 - r > 1. (40) 



Therefore the integrand in s(A) is a smooth function of r for r < n/2 and quickly decays to 
zero for r > n/2. The contribution to the integral in s(A) from the range of r <C n is small 
(^((r/n) 1 / 2 ). 



We note that term Fk in fl36| ) provides an "entropic" contribution to the sum in (|3g). It 
comes from the large number of states z G £k corresponding to large Hamming distances 
r from the state z', 1 < r < n/2. Each state contributes a small weight, G r oc ( n ) , and 
number of states for a given r is large, (cu^+i — k>fc) (™) -Pr 3> 1- Here {ujk+i — ujk) is an 
energy bin for the subspace and P r is the conditional density of states described in Sec. 
O. The size of the bin scales down exponentially with k (cf. (p7|)) and so does the entropic 
term F^. Below a certain cross-over value of k one has \Fk\ *C |/z',fe(A)|. In this case the 
dominant contribution to the sum ( p6|) comes from the states z with small r = -D(z, z') ~ 1. 
In particular for k = one can obtain 

/z',o(A) w Gx(A) ^ S lM<Vf) + 0(n- 3 ), (41) 

wG£o 



where the higher-order terms correspond to D(z', w) > 2. According to (39), |Gi(A)| ~ n~ 2 
and therefore |/ z ,o| is exponentially larger than the entropic term, |i*o| ~ ujq ~ d§2~ n . We 
note that, unlike the entropic term, / z / )0 strongly depends on z' due to the discreteness of 
the partition energy spectrum {oj^n ^ 1). E.g., depending on a state z', in this case there 
could be either one or none of the states w 6 £ in the sum ( pETD satisfying -D(z', w) = 1. 

B. Extended and localized eigenstates 

Based on the discussion above we look for solution of Eq. (^) in the following form: 

(f) z = v(Q z ) +u z , z<£ C , (42) 
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where we have explicitly separated out a part of the wavefunction v(Q z ) that depends on z 
only via the corresponding value of the partition residue. It satisfies the following equations: 



r$2~ n f°° 

[i-TG (x)c(Q)]v(n) = - tt + t dn'v(n')c(n')x(n',n,X), (43) 

A - aV 

x(n', n, A) = J2 ( n ) °rW P r (n'\n). (44) 



where $ is given in (g^) and function c(x) takes a set of discrete values (|6|). Using (p5|),(^2|) 
and fl43|) we obtain equations for u z 

M 

[1 - T G (A) E k ] U z = T ^ £ k' ^2 G D{z,z')W uw + TS ^ G D{z<w ) (A)0 W , Z G £fe. 

fc'=l z'e£ fc we£ 

(45) 



Decomposition (|42D is only applied to amplitudes <p z with z ^ C . The system of equations 
for the components v(Q) and u z is closed by adding Eq. (j35|) for the amplitudes W with 
w G Co (ground states of the final Hamiltonian Hp) and taking ( |42|) into account. We note 
that Eq. (|43"D for v (Q) is coupled to the rest of the equations only via the symmetric term $ 

$ = $ + $ + $ (46) 

/oo 
dxP(x)v(x)c(x), (47) 
-00 



M 



where distribution P(O) is given in 



Minimum gap estimate for ujm "C 

We will analyze the above system of equations ([42|)-(|47|) assuming that the cutoff fre- 
quency ujm satisfies Eq. fl29|) . This condition corresponds to the linear region in the plot of 
the cumulative density of states given in insert to the Fig. [II A| . According to Eqs. (^),(|19|). 
in this range the distribution functions -P(fi) ~ P n /2 an d P r (0'|0) ~ P r take nearly constant 
values and spectral function x(^', A) equals 

x{a,n,\)* /9 s(A) 2rnV (48) 
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where s(A) is given in (PT|). In this approximation, we can compute $ using equations for 
u z in ( ^5| ) and also the relations in (PE|), ( |5TD 

$ = -k(t/is(A))$ , k(x) = — — . (49) 

1 + x 

In the initial stage of the algorithm the amplitudes W of the "solution" states are small 
|$o| = 0(2~ n/2 ). According to ([49]), we also have |$| = 0(2~ n/2 ). Neglecting these terms 
and setting <3> « $, Eq. (^3|) gives a closed- form algebraic equation for A 

1 + 2 ^(a3W + S(A) ) = °' <50) 
Expanding in a small parameter /i < 1 (cf.(|29]),(|38|)), we obtain the eigenvalue 

A l (r) « a(r)F - 2r/i - 2(r/x)2 ln2 + 0(/i 3 ) (a » /,), (51) 

a 

that accurately tracks the adiabatic ground state energy, Ao(t), from r = 0, up until small 
vicinity of the avoided-crossing, r m r* (see below) where |$o| ~ 1- 

In the avoided-crossing region, branch Aq(t) intersects with another branch, Aq(t), that 
tracks Ao(t) in the interval of time following the avoided-crossing, r* < r < 1. This branch 
corresponds to $ <C $ , It can be obtained from simultaneous solution of equations for u z 
(fHf) and </> w that are approximately decoupled from Eq. (fy|) after $ is neglected. Keeping 
this term in ([15]) gives rise to repulsion between branches Ag (t) at r = r* that determines 
the minimum gap g m i n (see below). 

To proceed, we obtain the equation for <3> by adding equations for amplitudes W that 
correspond to different states w 6 £ and neglecting the coupling between these states 
separated by large Hamming distances, _D(w,w') ~ n/2. It can be shown using Eqs. fl35|) 
and (]4T|)-(f45|) that u z enters equation for $ through the term 

r 2 e Q ^2 ^ z / z ,o(A)m z , (52) 

which is is a self-energy term corresponding to elementary bit-flip processes with initial and 
final states belonging to the subspace £ (loop diagrams). 

To express u z in (|52|) through <j) w we solve Eq. ( f4"5j ) using order-by-order expansion in a 
small parameter n' 1 (cf. Eqs. (^)-(5T) and discussion there). In particular, one can show 



that to the leading order in n 1 the self-energy term (B2I) is determined by lowest-order loops 
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with two bit flips that begin and end at Cq. Then after some transformations, the equation 
for $o takes the form 

*° - - - ^ £ = ^02- + * W) . (53) 

Here a = 1 — r (cf. ( p4|) and $ is defined above. We now solve Eq. ([53]) jointly with (f43"D and 
obtain a closed- form equation for A. We give it below in the region of interest \r — 1/2 1 <C 1 

(A - A> (r)) (A - AS(r)) = -n 2 2-A 2 /4 (54) 
A « 4 /2 (1 + /ir* In 2 + 0(/i 2 )) , 

where the branch Aq(t) is given above and the branch Aq(t) satisfies Eq. (|53| ) with r.h.s. 
there set to zero, 

A f (r) « T£ - 1/2, |r-l/2|<l. (55) 
Avoided-crossing in (|54]) takes place at r = r* 

A l (r*) = Aj(r*), T^ + i-log,^. (56) 
The value of minimum gap between the two roots of ( |54]) equals 

ta = «A2- li / 2 . (57) 

where A is defined in (^7|). 

Based on the above analysis one can also estimate the matrix element \(4>i\H T \<f)o)\ T=T * ~ 
n. Then from Eq. (|25|) (see also discussion after Eq. (|28|) ) one can estimate the run-time of 
the quantum adiabatic algorithm 

T » rfo| f* Ql1 = O((7id ) _1 2 n ). (58) 

drain 

It follows from the above that eigenvalue branch Aq(t) corresponds to a state, 

|0 O }« £ «(n„)|z>, 

ze{o,i} n 

which is extended in the space of the bit configurations |z): according to ( f4"3"D it contains a 
large number (0(2™)) of exponentially small (0(2~™/ 2 )) individual amplitudes. This state 
originates at r = from the totally symmetric initial state 1^(0)) (0). In the small region 
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\t — r*| ~ g mm it is transformed into the state that corresponds to the eigenvalue branch 
Aq(t) and is localized in Hamming distances D(z,w) near the subspace w e £q containing 
the solution to NPP |0 O ) X] w e£ l w ^' Minimum gap at the avoided-crossing is determined 
by the overlap between the extended and localized states. 

At later times r > r* a similar picture applies to the avoided crossing of the extended- 
state energy Aq(t) with energies of localized states A^(t) corresponding to z 6 4 with 
1 < k <§C n (excited levels of the final Hamiltonian Hp (|28| ) ) . The existence of the extended 
eigenstate of H(t) whose properties do not depend on a particular instance of NPP follows 
directly from Eq. (|43|) that involves only a self-averaging quantity x(^'> ^> This quantity 
varies smoothly over the broad range of partition residues \Q'\, \Q\ < (E) and does not 
allow for the compression of the wave-packet v(Q z ) on the much smaller scale 0(2~ n ). This 
gives rise to an eigenstate with probability amplitude of individual states |z) that depends 
smoothly on energy in this range. 



2. Analysis of the general case 

The above picture of avoided-crossing remains qualitatively the same when the condition 
is relaxed (cf. insert in the Fig. ^). Away from the avoided-crossing point, r < r*, 
the ground state wavefunction v(Q z ) and energy Aq(t) are obtained directly from Eq. ([43]) 
with replacement $ ~ $ and Eq. ( pH| ) taken into account. Because the spectral function 
x(f2, A) changes only slightly on the scale E min = 0{n 1 / 2 2' n ) the wave packet v (f2 z ) |z) 
remains extended, \v{Q, x \ = C(2 _n/2 ), and therefore $ = 0(2~ n/2 ). 

Beyond the avoided-crossing point, r > r*, the ground state is localized near w and 
eigenvalue branch Aq(t) is obtained from Eq. (53) with r.h.s. set to zero (cf. Sec. |I V B 1| ) . 



The point r = r* is located at the intersection of the two branches A (r) ~ Aq(t) and the 
level repulsion is of the order of the overlap factor between the extended and localized states 



Ground-state wavefunction <p z at the avoided-crossing is shown in Fig. [| for modest value 
of n, but the separation into slowly- and rapidly-varying parts (|42]) is clearly seen. 

We did not perform a direct numerical study of the dependence of g min on n since we 
only simulated adiabatic eigenvalues for small instances of NPP. We argue, however, that 
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FIG. 2: Dots correspond to the plot of the ground state amplitude (z|0o) v $ partition residue 
\Q Z \ evaluated at the avoided crossing point r = r* (thin lines connecting the dots are for display 
purposes). Simulations are done for the randomly sampled instance of NPP with n = 10 and 
b = 20; the corresponding value of r* ~ 0.5. In simulations we relax the condition ( |29|) and the 
value of M in (27) is set automatically to be an integer closest to log 2 J2j=o a j ( c ^- (P^D)- Insert: 
Dotted curves are the plots of the two lowest eigenvalues of H (r) vs r for the same instance of 
NPP as in the main figure. Solid lines that start at r = correspond to A = (1 — r)n + k with 
k = 0, 1 (cf. d5l])). Solid lines that ends at r = 1 correspond to A = re^ with k = 0, 1 (cf. (|55|)), 



even for a fixed n the scaling of g m j n with n can be inferred from the shape of the cumulative 
density of states 

p X km 

77(A) = / dxJ2H^-x), k m = 2 n -l, (60) 
Jo k=0 

where A& = Ajt(r) are eigenvalues of H{r) (|4]). These eigenvalues are plotted in Fig. [|near 
the avoided-crossing r = r* where the spectrum of A& is quasi-continuous. The shape of the 
plot is well approximated by the square-root function: 

1/2 



X v ~ const + 



Via 



Vn 



0(2 n ). 



(61) 



It is clear that for 77 « 1 we have A^ ~ 2 ™/ 2 which corresponds to Eq. ([571). Note that 
this qualitative analysis is based on the assumption that the asymptotic properties of Aq for 
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FIG. 3: Dotted line is plot of A^ vs r] at the avoided-crossing r = r*. It is obtained from the 
numerical solution of the stationary Schrodinger equation for the same instance of NPP as in 
Fig. ^. Solid line is a square-root fit A = —6.3 + 0.35 7/ 1 / 2 (solid line is almost undistinguishable 
from the dotted line). 

large n can be inferred from the behavior of A^ for r\ 3> 1. 

V. SIMULATIONS OF TIME-DEPENDENT SCHRODINGER EQUATION 

We also study the complexity of the algorithm by numerical integration of the time- 
dependent Schrodinger equation with Hamiltonian H(t) and initial state \ip{0)) defined 



in Eqs. (p2|) ,(p3[),(p7D,(28). Here we relax the condition lum <^ (E) used above in the 
analytical treatment of the problem; in simulations the value of M is set automatically to 
be an integer closest to log 2 2^j=o a i ( c ^- We introduce a complexity metric for the 

algorithm, C(T) = (1 + T)d /p {T) where p {t) = J2 wECo l<Mt)| 2 . A typical plot of C(T) 
for an instance of the problem with n=15 numbers is shown in the insert of Fig. 4. At 
very small T the wavefunction is close to the symmetric initial state and the complexity is 
~ 2 n . The extremely sharp decrease in C(T) with T is due to the buildup of the population 
Po(T) in the ground level, S z = Eq, as quantum evolution approaches the adiabatic limit. At 
certain T = T min the function C{T) goes through the minimum: for T > T min the decrease 
in the number of trials do/po(T) does not compensate anymore for the overall increase in 
the runtime T for each trial. For a given problem instance the "minimum" complexity 
Cmin = C(T m i n ) is obtained via one dimensional minimization over T. The plot of the 
complexity C min for different values of n in Fig. 1 appears to indicate the exponential 
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scaling law, C m i n ~ 2°' 8n for not too small values of n > 11. 

1e+04 

1e+03 

Q 

min 

1e+02 



1e+01 

8 9 10 11 12 13 14 15 16 17 

n 

FIG. 4: Logarithmic plot of C m ; n vs n for randomly generated instances of NPP with 25-bit 
precision numbers. Vertical sets of points indicate results of different trials (~ 100 trials for each 
n, except n=17 with 10 trials). Median values of C m i n are shown with rectangles. Linear fit to the 
logarithmic plot of median values for n between 11 and 17 is shown by the line and gives lnC m i n ~ 
0.55n (C m i n ~ 2 8n ) . Very close result is obtained for the linear fit if all data points are used 
instead of the median values. Insert: plot of C(T) vs T for n=15, precision b=25 bits, do=22. 
Point 1 indicated with the arrow refers to the minimum Vctluc of complexity cit T — 2~min — 

22.67 

where the total population of a ground level po(^min) = 0.15. Point 2 refers to the value of T where 
Po (T) = 0.7. 

VI. DISCUSSION 

In conclusion, we have developed a general method for the analysis of avoided-crossing 
phenomenon in quantum spin-glass problems and used it to study the performance of the 
quantum adiabatic evolution algorithm on random instances of the Number Partitioning 
problem. This algorithm is viewed as a "quantum local search" with matrix elements of 
the Green function G r (r = 1, . . . , n — 1) giving the quantum amplitudes of the transitions 
with different number of spin flips r. Our approach is similar to the analysis of a quantum 
diffusion in a disordered medium with the model of disorder defined by the one- and two- 
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point distribution functions P(fl), P r . )Z (0 / ). 

We have shown that the conditional distribution of partition residues P r (Q'\Q) in the 
neighborhood of a given string formed by all possible r-bit flips depends on the value of the 
partition residue for that string but not on the string itself. This is a specific property of 
the random Number Partitioning problem. 

We used the above property to describe a quantum diffusion in the energy space 
(Eq. (|43|) ). This reduction in the dimensionality leads to the formation of the eigenstate 
which is extended in the energy space. Near the avoided-crossing the adiabatic ground state 
changes from extended to mostly localized near the solution to the optimization problem. 
Because the extended and localized state amplitudes are nearly orthogonal to each other the 
repulsion between the corresponding branches of eigenvalues (the minimum gap) is expo- 
nentially small, g m i n ~ n2~ n / 2 , and the run time of the algorithm scales exponentially with 
n. Analytical results are in qualitative agreement with numerical simulations of the time- 
dependent Schrodinger equation for small-to-moderate instances of the Number Partitioning 
problem (n < 17). 

One can show that the effect of quantum diffusion in reduced-dimensional space that 
leads to the formation of the extended state can also occur in other random NP-complete 
problems The method developed in this paper will be applied to study the performance 
of continuous-time quantum algorithms for different random combinatorial optimization 
problems. Also the present framework can be applied to the analysis of quantum annealing 
algorithms for combinatorial optimization problems []32| , j33f . This is a classical algorithm 
that is conceptually very close to the quantum adiabatic evolution algorithm considered 
above |54j]. The former uses the Quantum Monte Carlo method to simulate on classical 
computers a partition function and ground-state energy of a quantum system with slowly 
varying Hamiltonian that merges at the final moment with the problem Hamiltonian of a 
given classical optimization problem. Among other possible applications of our method is 
the analysis of tunneling phenomenon in the low-temperature dynamics of random magnets. 

We note that the specific property of the Number Partitioning problem (that distinguishes 
it from the other NP-complete problems) is a very weak dependence of P r (Q'\Q) on Q for not 
too large values of <C ^Jr(n — r) that takes place for all values of r e [1, n — 1]. This 
rapid fall-off of correlations during the local search (both classical and quantum) is a reason 
that the exponential complexity of optimization algorithms for the Number Partitioning 
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problem can be seen already for the relatively small values of n < 15 (cf. FigJJ). 

Finally, our analysis of sub-harmonic resonances in the Fourier transform 7(s) of the 
distribution function P(fi) suggests a possible connection between NPP and the integer 
factorization problem. If, for a given set of Oj's, there is a number q that satisfies the 
condition ( |A5| ) then dividing all numbers aj by q we obtain a new instance of NPP with 
numbers kj = ctj/q that will be completely equivalent to the old one. It is important that 
the precision of the numbers kj is restricted by b — log 2 q. If the value of q is sufficiently 
large, log 2 q ^> b — n, then k^s correspond to a low precision instance of NPP, i.e. to 
the computationally easy phase mentioned in Sec. y. This is exactly the case when sub- 
harmonic resonances become substantial. One can fix the parameter £ = b/n ^> 1 in a 
high-precision (computationally hard) case and compute, for randomly generated instances 
{cij} an approximate greatest common divider, i.e. a largest number q that satisfies (|X3|). The 
distribution of these numbers determines a fraction of high-precision instances of NPP (out 
of all possible 2 nb problem instances) that really belong to a low-precision (computationally 
easy) "phase". 

Advance knowledge of this information would be of importance if one is using NPP 
for encryption purposes flI7| , especially because NPP is otherwise a very difficult problem 
for both quantum and classical computers It is not obvious at this stage what the 

asymptotic form of this distribution will be in the limit of large n (cf. Fig. |AT). 

We are not aware of any classical algorithm that could verify if such a number q exists 
for a given set of aj in a time polynomial in both n and b. However, on a quantum computer 
one can apply a Shor algorithm to test in polynomial time if strong sub-harmonic resonances 
exist. This question is deferred to a future study. 
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APPENDIX A: SUB-HARMONIC RESONANCES 




We note that function J(s) in (]5|) can also have additional sharp resonances in the range 
< \s\ < 2 b . To understand their origin we consider first a particular case when rational 
6-bit numbers 01,0,2, • • . , a n all have a number q > 2~ b as a "common divisor", i.e., there 
exist integers ki, A; 2 , . . . , k n such that 

— — — — — — — q (Al) 

h k 2 k n 

In this case additional resonances of I(s) occur at the multiples of n/q. Assume now that 
q is no longer an exact divisor of numbers aj but all the residues of the divisions a,j/q 
are sufficiently small. Then contributions from the additional resonances at s ~ mir/q 
(m = 1,2,...) to the integral in (|5]) can be estimated as follows (for simplicity we give a 
result for the case E <C n 1 ^ 2 ): 

on 

P(0) -> y= > ] C [ ) (-l) mp (A2) 

P = 

Here [x] and {x} denote integer and fractional parts of a number x, respectively. If the total 
"dephasing" factor e^ 1 ^ ~ 1, then contribution ( |A2| ) cannot be neglected in the steepest- 
descent analysis of (^) (in general, on should keep contributions from all divisors q with 
small dephasing factors e -7 ^). 

We note that the window function ( (^^j ~ 1 for g ^> 2~ n and it decays to zero 
at smaller values of q. We studied numerically the greatest root g max of the the algebraic 
equation 

i(q) = ic (A3) 

for a fixed value of 7 C < 1. For the sets of random 6-bit numbers aj the dependence of 
the mean value of q m { n on the problem size n < b is shown in Fig. |5|. For n <ti b we have 
exponential decrease of g max with n and for larger values of n < b the value of q m \ n steeply 
drops to 1. According to the discussion above, in order to neglect the saddle-points with 
s > in (|5p (additional resonances) the value of g m ; n should satisfy the following condition 
in the asymptotic limit b — > 00: 

< max [2~ n , 2~ b ] , 1< n < b, (A4) 
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with 7 C fixed at some small constant value. Because the precision b that we used in the 
simulations was not very high (limited by machine precision) it is not possible to obtain the 
asymptotic form of the dependence of q m [ n on n in the range given in (|A4[ ). Neither we can 
describe the shape of the plot in Fig. |5] analytically in that range. However, it appears from 
the figure that the condition ( |A4| ) is satisfied for sufficiently large n. 




n 



FIG. 5: Log-Log plots of the mean value of the largest root of Eq. (A3) g m ; n vs n. Three sets 
of data points are plotted. Each set of points represents averaging over 25 randomly generated 
instances of NPP. Precision of the random numbers aj is 30 bits and the value of 7 C = 0.5. Dashed 
line corresponds to the plot of const x 2 _n vs n. Insert: Variance of the log 2 q vs n based on 25 
sample points for each n. Distribution of q m i n values become very broad when the mean drops to 

<?min ~ 1- 



APPENDIX B: PROPERTIES OF THE CONDITIONAL DISTRIBUTION OF 
SIGNED RESIDUES IN NPP 

We perform the summation over the spin configurations in Eq. (|T7|) with Eq. (|l~5D taken 
into account. Similar to the derivation of Eq. @ we use integral representation for delta 
function and obtain 

rw^A f°° r dsds' fAQs\ /AfiV\ ^ TT/ , N 

Mv-U LL^ ( {—) c {—)^ u ^ s) ' (B1) 

Uj(s,s') = Y[coB(aj(8- 8 ')) x ^005(^(3 + x e i(sn+s ' n,) . 
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Here the sum is over all possible subsets J = {ji,j2, ■ ■ ■ ,jr} of length r obtained from the 
set of integers j — 1, 2, . . . , n. Window function Q(x) is defined in ([5]). After the change of 
variables 

x' = s + s\ x = s- s', (B2) 
we obtain from ( |B1| ) that Uj(s, s') factorizes into a product of two terms 

U 3 (s,s') = Vj(x)Vj(x') 



Vj W = exp ( I || cos{aj x), Vj{x ) = exp I I || cos(aj x ). 



(B3) 



In what follows we will analyze several limiting cases. 



r, n — r ^> 1: 



In this case both functions Vj(x) and Vj(x') are very steep and similar to the analysis 
in Sec. [II A| integrals in (|B1|) can be evaluated by the steepest descent method. With the 
appropriate choice of the coarse-graining windows AQ, AQ' in ([Bl]) (see below) contribution 
to the integrals comes from the vicinity of the point (x — 0, x' — 0). Near this point we use 



J J cos(aj x) ~ exp 



| J cos(a.,- x) « exp 



(n - r)(x / crj' 



(B4) 



where 



n — r 



jeJ j£j 
Since each sum here contains a large number of terms we obtain for i.i.d. random numbers 

ai,...,a n (cf. (I)) 



(<7j) 2 «<t 2 (O)+0 - 



'aj) 2 ^(O) + 



1 



(B5) 



r / \n — r / 

where a 2 = (a 2 ) is given in (||). Using Eqs. (p3|)-(p5|) and replacing the window functions 
in (|B1| ) by unity, we compute the Gaussian integrals in ( pi]) and obtain 

1 



47rcr 2 (0) \/r(n — r] 



exp 



?<r 2 (0) 



i f(n-n') 2 (fi + ^') 2vl 



r 



n — r 



(B6) 



The size of the coarse-graining windows in ([Bl]) is chosen to satisfy the conditions 

T n {™\ < Aft Afi' < y/r(n - r) 
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From Eq. (|B6| ) and Eq. (|6]) one can directly obtain the conditional distribution function 



P r (Q'\Q) = -===e W 



2na 2 (q) 



(B7) 



n — 1: 



For r = 1 function Vj(x') contains a product of n — 1 terms and is very steep. The 
corresponding integral over x' in ( |ET| ) should be taken by the steepest descent method. 
However Vj(x) simply oscillates at frequencies (Q — Q')/2 ± a, and the integral over x in 
( |B1[ ) should be evaluated using the corresponding oscillating factors. In the opposite case 
r = n — 1, function Vj(sc) is very steep and the integral over x in ( |B1| ) should be taken 
by steepest descent. But the integral over x' there should be evaluated using Vj(x') that 
oscillates at the frequencies, (Q + Q')/2 ± dj. Finally, one can obtain using i.i.d. numbers 
aj's in [0, 1] interval : 

P r (Q'\Q) = i[e(fi=Ffi' + 2) -e(n=Ffi' -2)] + f^j , (r = l, n-1). (B8) 

The minus (plus) sign in (|B8| ) corresponds to r = 1 (r = n — 1). Similarly one can obtain 
the result for any fixed value of r or n — r (that does not scale with n). For < 1 

M) is reduced to (TO 



Numerical simulations of conditional distribution P r z (fi') 



We compute the following integrated quantity: 



1 



Q=2j o d VPrM, (B9) 

for different values of r, Q' and different strings z with i? z 1. Numerical results are 
compared in the insert to Fig.^Jwith theoretical result below obtained using P r (Q'\Q) from 
Eq. Q 

-/ eft ? P r (t ? |0)=erf . (BIO) 

Theoretical and numerical curves nearly coincide with each other. To accurately compare 
the normalization factor in (|B7|) (see also (0)) we compare the theoretical results with 
numerical values of P r , z (0) for different r and strings z corresponding to E z 1. The 
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results are plotted in Fig. [1| 



[1] P.W. Shor, in Proceedings of the 35th Annual Symposium on the Foundations of Computer 
Science, ed. by S. Goldwasser (IEEE Computer Society Press, Los Alamitos, CA, 1994), p. 124; 
SIAM J. Comput. 26, p.1484 (1997). 



[2] 
[3] 

K 

[5] 
[6] 

[7] 
[8] 
[9 



[10 

[11 
[12 

[13 
[14 

[15 

[16 

[17; 



Y. Fu and P.W. Anderson, J. Phys. A: Math. Gen. 19, 1605-1620 (1986). 
M. Mezard, G. Parizi, and M. Virasoro, Spin glass theory and beyond (World Scientific, Sin- 
gapore, 1987). 

M.R. Garey and D.S. Johnson, Computers and Intractability. A Guide to the Theory of NP- 
Completeness (W.H. Freeman, New York, 1997) 

E. Farhi, J. Goldstone, S. Gutmann, and M. Sipser, |arXiv:quant-ph/0001106 . 



E. Farhi, J. Goldstone, S. Gutmann, J. Lapan, A. Lundgren, and D. Preda, Science 292, 472 
(2001). 

E. Farhi, J. Goldstone, and S. Gutmann, arXiv:quant-ph/0007071. 



A. M. Childs, E. Farhi, J. Goldstone, and S. Gutmann, |arXiv:quant-ph/0012104 



T.Hogg, 'Adiabatic Quantum Computing for Random Satisfiability Problems" , arXiv: quant- 



ph/020605£ 



W. Van Dam, M. Mosca, U. Vazirani, "How Powerful is adiabatic Quantum Computation?", 
|arXiv:quant-ph/020600~3 . 



C. Bennett, E. Bernstein, G. Brassard, and U. Vazirani, "Sterngths and weaknesses of quantum 
computing", SIAM Journal of Computing, 26, pp. 1510-1523 (1997); |arXiv:quant-ph/9701001 



W. Van Dam, M. Mosca, U. Vazirani, "How Powerful is Idiabatic Quantum Computation?", 
FOCS 2001. 

E. Farhi, J. Goldstone, S. Gutmann, |arXiv:quant-ph/0201031 . 



A. M. Childs, E. Deotto, E. Farhi, J. Goldstone, S. Gutmann, A. J. Landhal, "Quantum 
search by measurment" , |arXiv:quant-ph/0204013 . 



A. Boulatov, V. Smelyanskiy, "Total suppression of a large spin tunneling barrier in quantum 
adiabatic computation", arXiv:quant-ph/0208189| . 



Li-Hui Tsai, SIAM J. Comput., 21(1) p.59-64 (1992). 

A. Shamir, Proc. of 11th Annual ACM Symposium on Theory of Computing, p. 118 (1979). 

30 



[18] N. Karmakar, R.M. Karp, G.S. Lueker, and A.M. Odlyzko, J. App. Prob. 23, p. 626 (1986). 

[19] D.S. Johnson, et al., Operations Research 39, p.378 (1991). 

[20] F.F. Ferreira and J. F. Fontanari, J. Phys. A 31, p. 3417 (1998). 

[21] LP Gent and T. Walsh, in Proc. of the ECAI-96, ed. by W. Wahlster (John- Wiley & Sons, 

New York, 1996), pp. 170-174. 
[22] R.E. Korf, Artificial Intelligence 106, 181 (1998). 
[23] S. Mertens, Phys. Rev.Lett. 81, 4281-4284 (1998). 

[24] C. Borgs, J. T. Chayes and B. Pittel, Proc. of the 2001 ACM Symposium on the Theory of 

Computing, pp. 330-336 (2001). 
[25] C. Borgs, J. T. Chayes and B. Pittel, Random Structures and Algorithms, v. 19, pp. 247-288 

(2001) 

[26] S. Mertens, "A complete anytime algorithm for balanced partitioning" , arXiv:abs/cs.DS/9903011. 
[27] P. Cheeseman, B. Kanefsky and W. M. Taylor, Proc. of the International Joint conference on 

Artificial Intelligence, v. 1, pp. 331-337 (1991). 
[28] Artif. Intel. 81 (1-2) (1996), special issue on Topic, ed. by T. Hogg, B.A. Huberman, and C. 

Williams. 

[29] S. Mertens, Phys. Rev.Lett. 84, 1347-1350 (2000). 
[30] S. Lloyd, Science 273, 1073 (1996). 

[31] V. Smelyanskiy, D. Shumow, U. Toussaint, in preparation. 
[32] T. Kadowaki and N. Nishimori, Phys. Rev. E 58, p. 5355 (1998). 
[33] G. Santoro, R. Martonak, E. Tosatti, R. Car, Science 295, p. 2427 (2002). 
[34] T. Kadowaki, "Study of Optimization problems by quantum annealing", PhD Thesis, 
arXiv:quant-ph/020520 (2002). 



31 



